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Abstract 

Texture segmentation constitutes a standard image processing task, crucial for 
many applications. The present contribution focuses on the particular subset of scale- 
free textures and its originality resides in the combination of three kc!y ingredients: 
First, texture characterization relies on the concept of local regularity ; Second, esti¬ 
mation of local regularity is based on n(!w multiscale quantities referred to as wavcdet 
leaders : Third, segmentation from local rc!gularity faces a fundamental bias variance 
trade-off: In nature, local regularity estimation shows high variability that impairs the 
detection of changes, while a posteriori smoothing of regularity estimates precludes 
from locating correctly changes. Instead, the present contribution proposes several 
variational problem formulations based on total variation and proximal rc!Solutions 
that effectively circumvent this trade-off. Estimation and S(;gmentation performance 
for the proposed procedures arc; quantified and compared on synthetic as well as on 
real-world textures. 


1 Introduction 


Fractal based texture characterization. Texture characterization and segmentation 
are long-standing problems in Image Processing that received significant research efforts in 
the past and still attract considerable attention. In essence, texture consists of a perceptual 
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attribute. It has thus no unique formal definition and has been envisaged using several 
different mathematical models, mostly relying on the definition and classification of either 
geometrical or statistical features, or primitives (cf. e.g., and references therein for 

reviews). Texture analysis can be performed using either parametric models (ARMA 
Markov Wold j^) or non parametric approaches (e.g., time-frequency and Gabor dis¬ 
tributions iHD- Amongst this later class, multiscale representations (wavelets, contourlet, 
...) have repeatedly been reported as central in the last two decades (cf. e.g., IHnl). 
They notably showed significant relevance for the large class of scale-free textures, often 
well accounted for by the celebrated fractional Brownian motion (fBm) model, on which the 
present contribution focuses. Examples of such textures are illustrated in Fig. [^a). While 
scale free-like texture segmentation has been mainly conducted by exploiting the statistical 


distribution of the pixel amplitudes 14 , the present paper investigates the relevance of 
local regularity-based analysis. 


Deeply tied to multiscale analysis, the fractal and multifractal paradigms |15[[T6| have 
been intensively used for scale-free texture characterization (cf. e.g., [T7[[l8| ). Scale-free 
textures are essentially measured jointly at several scales, or resolutions. From the evolu¬ 
tion across scales of such multiscale measures, semi-parametrically modeled as power-laws 
and thus insensitive to texture resolution, scaling exponents are then extracted and used as 
characterizing features. Such fractal features have been extensively used for texture char¬ 
acterization, notably for biomedical diagnosis (e.g., [T& 
investigations [^[^. 


22 ), satellite imagery 1231 or art 


Often, in applications, it is assumed a priori that fractal properties are homogenous across 
the entire piece of texture to characterize, thus permitting reliable estimates of the fractal 
features (scaling exponents). However, in numerous situations, images actually consist of 
several pieces, each with different textures, and hence different fractal properties. Analysis 
becomes thus far more complicated as it requires to segment the texture into pieces (with 
unknown boundaries to be estimated) within which fractal properties can be considered 
homogeneous (i.e., fractal attributes are comstant), yet unknown. 


Local regularity and Holder exponent. To address such situations, the present contri¬ 
bution elaborates on the fractal paradigm, by a local formulation relying on the notion of 
pointwise regularity. The local regularity of a function (or sample path of a random field) 
at a given location r € is most commonly quantified by the Holder exponent h{x) (^ . 
It is defined as the scaling exponent extracted from the a priori assumed power law depen¬ 
dence of the coefficients of a multiscale representation X(a, r) (e.g., the modulus of wavelet 
coefficients 11 ), across the scales a, in the limit of fine scales: X{a.x) ~ T]{x)a^^^ when 
a ^ 0. Theoretically, the collection of Holder exponents h{x) for all x G yields access to 
the multifractal properties (and spectrum) of the texture, and thus consists of a multiscale 


higher-order statistics feature characterizing the texture 26 27 


Wavelet coefficients are the most popular multiscale representation used to perform fractal 
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analysis. Yet, it has recently been shown that, in the context of multifractal analysis and 
thus of local regularity estimation, wavelet coefficients are significantly outperformed by 
wavelet leaders^ consisting of local suprema of wavelet coefficients [26| - [^ . While the meth¬ 
ods proposed in the present contribution could be used with any raultiscale representation, 
reported and discussed results are thus explicitly obtained using wavelet leaders. 

For each location h (^) is classically estimated via a linear regression in log X (a, versus 
log a coordinates and is thus naturally framed into a classical bias-variance trade-off: Point- 
wise estimates (relying only on the X(a,^) at sole location show very large variances. 
Conversely, averaging in space-windows either the X{a,x) or directly the estimates of h{x) 
results in significant bias. In either case, the reliable and accurate detection and localization 
of actual changes in h is precluded. This explains why local regularity remains so far barely 
used for signal or texture characterization and segmentation (see, a contrario, |29{-|^). 


Goals, contributions and outline. In the present contribution, elaborating on prelini- 

the overall goal is to enable the actual use of local regularity for 


inary attempts 129 32 


performing the segmentation of scale-free textures into areas with piece-wise constant fractal 
characterization. This strategy has the great advantage of being fully nonparametric, since 
it does not require any explicit (statistical) modeling of the textures to be segmented. To 
that end. it aims to marry wavelet-leader based estimation of local regularity (with no a 
priori smoothing) to segmentation procedures formulated as variational problems and con¬ 
structed on total variation (TV) optimization procedures and proximal based resolutions. 
Wavelet-leader based estimation of Holder exponents (as defined and detailed in Section 
is first applied locally throughout the entire image. Then, to go beyond a naive a posteriori 


smoothing and threshold-based segmentation procedure (as described in Section 3.11, three 


different TV-based .segmentation procedures are theoretically devised: i) Local estimates of 
h are subjected a posteriori to a TV based denoising procedure, followed by a thresholding 


step for segmentation (cf. Section 3.2); ii) Local estimation of h is a priori embodied into 
the TV based procedure aiming to favor piecewise constant estimates, thresholding for seg¬ 


mentation is performed a posteriori (cf. Section 3.3): hi) Local estimates of h are subjected 
a posteriori to a TV based segmentation procedure, avoiding the denoising/thresholding 
steps (cf. Section 3.41. This later procedure is inspired by the convex relaxation of the 


Potts model detailed in 

In Sectionj^ the performance of the proposed TV segmentation procedures are compared 
and discussed for samples of synthetic textures with piece-wise constant h, obtained from 
a multifractional model, whose definition is customized by ourselves to achieve realistic 


textures (described in Section 2.4). Several different geometries and different sets of values 


for the regularity h of the different segments are considered. Furthermore, the impact of the 
TV-regularization parameter is investigated. These procedures are also compared against 
two texture segmentation procedures chosen because they are considered state-of-the art 
in the dedicated literature |^[^. The proposed approaches are further shown at work 
on samples chosen within well accepted references texture databases, such as the Berkeley 


3 









Segmentation Dataset. Synthesis and analysis procedures will be made publicly available 
at the time of publication. Note that the present contribution significantly differs from the 
preliminary works in that several methods involving TV are designed, and these 

methods are unified and compared on synthetic images and real-world textures. 


2 Holder exponent and wavelet leaders 


2.1 Holder exponent and Local regularity: Theory 


Let / = {f{x)) w'ith X = (xi, X 2 ) S fl, denote the bounded 2D function (image) to analyze. 
The local regularity around location x^ is quantified using the so-called Holder exponent 
/ 2 .(xq), formally defined as the largest a > 0 , such that there exist a constant x > 0 f^d 
a polynomial of degree lower than q, such that ||/{x) — 'Pin(x)|| < xll^ “ ^ 

neighborhood x of Xq, where || • || denotes the Euclidean norm 1^. When /i(xq) is close to 
0, the image is locally highly irregular and close to discontinuous. Conversely, large values 
of correspond to locations where the field is locally smooth. An example of a texture 

with piece-wise constant function h{x) is illustrated in Fig. (a). 


Though the theoretical definition of the Holder exponent above can be used for the 
mathematical study of the regularity properties of fields, it is also well-known that it turns 
extremely uneasy for practical purposes and for the actual computation of h from real world 
data. Instead, multiscale representations provide natural alternate ways for the practical 
quantification of h 27 ^ 30 . It is however nowadays well documented that the sole 
wavelet coefficients do not permit to accurately estimate local regularity 26 Early 


contributions on the subject proposed to relate local regularity to the skeleton of the con¬ 
tinuous w'avelet transform 11 23 . which however does not permit to achieve estimation of 
h at each location x as targeted here. Instead, it has been recently shown that an efficient 
estimation of h can be conducted using wavelet leaders 26 They are defined as local 
suprema of the coefficients of the discrete wavelet transform (DWT) and thus inherit their 
computational efficiency. 


2.2 Holder exponent and wavelet leaders 

Wavelet coefficients. Let 0 and ^ denote respectively the scaling function and mother 
wavelet, defining a ID multiresolution analysis [^. The mother wavelet 7 /; is further char¬ 
acterized by an integer Njp > 1, referred to as the number of vanishing moments and defined 
as: VA: = 0.1, ..., = 0 and f \t\^’fip{t)dt 7 ^ 0. From these univariate 
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functions, 2D wavelets are defined as 


j/;(2)(^) = {I){xi)tl’{x2), 1p^''^\x) = ■tp{xi)i}{x2). 

The 2D-DWT (L^-normalized) coefficients of the image / axe defined as 


( 1 ) 


( 2 ) 


where = 2 ^x — k),j e N*,fc € N^,m = 0,1,2,3} is the collection of 

dilated (to scales a = 2^) and translated (to locations x = 2^k. k = {ki,k 2 )) templates of 
[x). Interested readers may refer to [TIJ for further details. 

Wavelet leaders. The wavelet leader Lf{j,k), at scale 2^ and location x = 2^k, is defined 
as the local supremum of all wavelet coefficients taken across all finer scales 

2^' < 2Y within a spatial neighborhood 26 -^ 


L^^\j,k)= sup 

m={l,2,3} 


\2^-fYj^\f,k'] 


(3) 


y C3A.j 




(4) 


where 

=[fc 2 ^(fc + l) 2 J), 

3Aj,fc = UpG{-i,o,i}2 

The additional positive real parameter 7 can be tuned to ensure minimal regularity condi¬ 
tions for the case where the image / to analyze can not be modeled as a strictly bounded 
2D-function. It is set to 7 = 1 for the present contribution and not further discussed. In¬ 
terested readers are referred to [27l|28| for the details on the role and impact of parameter 


that the Holder exponent h{x) at location x is measured by 

( 5 ) 


It has been proven in I 
wavelet leaders as 

L^j^j^k) — r){x)2^^^^~'~'^^ when 2 -^ —>■ 0 
for x € Xj^ki provided that is strictly larger than h{x). 


Note that ^ implies that h can be estimated by linear regressions as the slope of log 
versus j, cf., ([^-(j^ below (similarly, log 7/ could be obtained as the intercept, yet does not 
bear any information on local regularity and will hence not be further considered here). 


2.3 Holder exponent estimation 


In the present contribution, the estimation of the Holder exponent is only performed using 


wavelet leaders L^j\j,k), as preliminary contributions 


27 


29 
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report that w'avelet leader 
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based estimation outperforms those based on other rnnltiresolution quantities. To indicate 
that the estimation of h and the segmentation procedures proposed in Section below 
could be applied using any other multiresolution quantity, e.g., the modulus of the 2D- 
DWT coefficients, a generic notation X{j,k) is used instead of the specific Yet, 

all numerical results in Section are obtained with wavelet leaders, X{j,k) = 

In the discrete setting, Hblder exponents are estimated for the locations x = 2k associated 
with the finest scale j = 1 , and we make use of the notation when we are dealing with 
matrices, rather than with matrix elements, e.g., 

K = {l,...,Ari} X {L...,W2}. 

As a preparatory step, the multiresolution quantity X is up-sampled in order to obtain as 
many coefficients at scales j > 1 as at the finest scale j = 1 

X(j, {kuk^)) = XO', {[2-4-11, r2“42l)) (6) 

with 1 < ki < Ni. 1 < ^2 ^ ^ 2 - With a slight abuse of notation, we will continue writing 
X for the upsampled multiresolution coefficients X. 

The relation (|^ can obviously be rewritten as 

log2 k) ^ jh{k) + log2 Tlik) (7) 

which naturally leads to the use of (weighted) linear regressions across scales for the esti¬ 
mation of h{x) 

^4) = X^^0'4)log2Y0',fc). (8) 

j 

Combining Q and (j^ above shows that, for each location k, the weights w{j.k) must 
satisfy the following constraints to ensure unbiased estimation 

and '^jw{j,k) = l. (9) 

j 3 

Though unusual, let us note that the weights w{j,k) can in principle depend on location k. 
Tliis will be used in the segmentation procedure defined in Section |.3.3| An illustration of 
such unbiased estimates, with a priori cho.sen iu{j,k) that do not depend on location k, is 
presented in Fig. (c). 


2.4 Piece-wise constant local regularity synthetic processes 

To illustrate the behavior of the segmentation procedures proposed in Section as well 
as to quantify and compare segmentation performances from the different procedures, use 
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is made of realizations of synthetic random fields with known and controlled piece-wise 
constant local regularity. These are constructed as 2D nmltifractional Brownian fields 
which are among the most widely used models for mildly evolving local regularity. Their 
definition has been slightly modified here to ensure the realistic requirement of homogeneous 
local variance across the entire image 



fix) = Cix) 




n 1 ^ 


Mx)+i 


-dWiO 


( 10 ) 


where is 2D Gaussian white noise and h{x) denotes the prescribed Holder exponent 

function. The normalizing factor C{x) ensures that the local variance of / does not depend 
on the location x. The details of the models do not matter much for the present work, as we 
are only targeting the control of local regularity of the field f{x). In the present contribu¬ 
tion, the function h{x) is chosen as piece-wise constant. Numerical procedures permitting 
the actual synthesis of such fields have been designed by ourselves. A typical sample field 
of such a process is shown in Fig. (a). 


3 Local regularity based image segmentation 


In this section, we will detail the proposed procedures for segmentation in homogeneous 
pieces of textures with constant local regularity. 


3.1 Local smoothing 


A straightforward attempt for obtaining labels consists in thresholding the histogram of 
pointwise estimates obtained with (j^ and (j^. However, such a solution yields estimates 
with prohibitively large variance. This prevents the identification of modes in the histogram 
that would correspond to the different zones of constant local regularity, see Fig. (c) for 
an illustration. The variance can be reduced a posteriori by a local spatial smoothing with 
a convolution filter g\ 

~ - 
h = g * h 

where h denotes the estimate described in Section 


2.3 For instance, g can model a local 
spatial average. In this work, we will consider a Gaussian smoothing parametrized by its 
standard deviation tj. Note that particular cases of this smoothing can be expressed as a 
variational formulation 


h =arg min 



Yjx{j,k)log2X{j,k) -h(k)'^ +'^||r4l 


itllF 


( 11 ) 


3=n 


7 





where || • ||f denotes the Frobenius norm and the transform F and the regularization param¬ 
eter A > 0 are related to the convolution kernel g. In particular, when A = 0, the smoothed 
-~s . “ 

solution h reduces to its non-smoothed counterpart h. An example of such an estimate 
with a = 10 is given in Fig. (d). Clearly, the variance of h is smaller than the variance of 
h (see Fig. (c) and (d)). Yet, it remains hard to identify separate modes in the histogram. 
What is more, the local averaging introduces bias at the edges of areas with constant h 
values, hence prevents any accurate localization of the regularity changes in the image. 


3.2 TV denoising 


To overcome these difficulties, we propose the use of total variation (TV) based optimization 
procedures, naturally favoring sharp edges, instead of local spatial averages. It is well known 
that the bounded total variation space, that is the space of functions with bounded ^i-norm 
of the gradient, allows to remove undesirable oscillations while it preserves sharp featmes. 
Rudin et al. 39 formalize this recovery problem as a variational approach involving a 


non-smooth functional referred to as total variation. This has been widely used in image 


processing for image quality enhancement 39 40 although it has been noted that it is not 


always well suited for restoration purposes due to the piece-wise constant nature of the 
restored images. In the present context of detection of local regularity changes, precisely 
such a piece-wise constant behavior of the solution is desired. 


The corresponding minimization problem reads 

S'^^=argmin| i V ( fc) log2 X(j, - h{k)) + XTV{h) 


( 12 ) 


where 


ke/c j=ji 


Ni—l N2—I 


Tv(o = E E ^{imih,h)y + {{D,k)(h,h)y, 

^cj=l A ;2 = 1 

(Z)i/i)(fci, ^ 2 ) = h{ki -|- 1, A:2 -|- 1) — h(^ki + 1, ^ 2 ), 

(Z)2/l)(^i: ^ 2 ) = h{ki -|- 1, A:2 -|- 1) — h(^ki, k2 1). 

The weights are a priori fixed, are independent of k. w{j,k) = w{j). and satisfy the con¬ 
straints (j^. Here, A > 0 models the regularization parameter that tunes the piece-wise 
constant behavior of the solution. Several techniques have been proposed in the literature 
to solve (R^, see for instance (iOlliTj. A forward-backward algorithm (4^ applied on the 


dual formulation of (121 is used here. 


'TV 


An example of h is displayed in Fig. (e) for A = 6. It is piecewise constant and 
provides a satisfactory estimate of the true regularity displayed in Fig. (b). Furthermore, 
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Figure 1: Illustration of the local regularity estimation from a imiltifractional data. Top 
row: (a) image /; (b) the original local regularity h from which / is generated, the area in 


black (resp. white) corresponds to a local regularity of 0.5 (resp. 0,7); (c 


3.1 


an estimation 
e) estimation 


based on wavelet leaders; (d) smooth estimate of h described in Section 
using the method described in Section [3.2[ (f) estimation based on the iiiethod described in 
Section 3.3 Bottom row: histograms corresponding to top row. 


-TV 

the histogram of the estimates h is pronouncedly peaked. Thresholding hence enables 
labeling for a preset number Q of classes. 


3.3 Joint estimation of regression weights and local regularity 

-TV 

The solution h described in the previous section is a two-step procedure that addresses 
the bias-variance trade-off difficulty by: (i) computing unbiased estimates of the Holder 
exponents h from X using (j^ with fixed pre-defined weights w, and (ii) extracting areas with 

constant Holder exponent based on these estimates h using a variational procedure relying 
on TV. In this section, we propose a one-step procedure that directly yields piece-wise 
constant local regularity estimates h from the multiresolution quantities X. The originality 
of this approach resides, on one hand, in the use of a criterion that involves directly X , 

instead of the intermediary local estimates h, and a TV-regularization; on other hand, in 
the fact that the weights w are jointly and simultaneously estimated instead of being fixed 
a priori. To this end, the weights w are subjected to a penalty for deviation from the hard 
constraints (j^. 

Problem formulation. The estimation (j^ underlies a linear inverse problem in which 
the estimate h of h needs to be recovered from the logarithm of the multiresolution quantities 
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X of the image /. This inverse problem resembles a denoising problem, yet including the 
additional challenge that a part of the observations (the regression weights r^) are unknown 
and governed by the constraints (j^, yielding the following convex minimization problem: 


I .?2 k 

(S ,w) =argmin< V ( 'yw{j,k)\og^X{j,k) - 




.k^K. j=ji 


+ ATV(^) + dc^{’^t^{-,k)) + r]2 (13) 


k€lC 


k€lC 


where J = J 2 “ + !■ The functions dci <^C 2 distances to the convex sets 


J 2 


Cl = {(a;(ji),...,u;(j2)} e | = 0}, 


3=Jl 

h 


C 2 = {(^(ji),• • e I = 1 } 


(14) 


(15) 


3=J1 


and are defined as 

(Vv e Ri)(Vi = 1, 2) dc,iv) = II 2 - PcAii)h 

where Pciin) = argmin^^gQ lln ~ fienotes the projection onto the convex set Ci- Note 

that the distances dci and dc 2 provide the possibility to relax the hyperplane constraints Ci 

and C 2 : The choice r/i = 772 = imposes (j^ as hard constraints (i.e., the intermediary 

quantities k) log 2 k) are unbiased estimates of h{k)), while for 771,772 < + 00 , 

a violation of the constraints ® is possible but penalized. This adds a degree of freedom 

^ ^ -TV _ __ 

as compared to the standard estimation procedure ^8^ and the solution h in (12). Fm- 


thermore, note that the joint estimation of h and w enables the use of spatially varying 
weights w, which is otherwise impractical for a priori fixed weights (here, the weights are 
tuned automatically). 


Proposed algorithm. The minimization problem (131 is convex but non-smooth. In the 


recent literature dedicated to non-smooth convex optimization, several efficient algorithms 
have been proposed, see, e.g., [43] - [45| . Due to the gradient Lipschitz data fidelity term and 
the presence of several regularization terms (TV regularization, distances to convex sets), 


one suited algorithm is given in [41 46 48], referred to as FBPD (for Forward-Backward 


Primal-Dual). This algorithm is tailored here to the problem (131, ensuring convergence 
of a sequence (h^^,u^^^)£eN to a solution of (13). The corresponding iterations are given in 
Algorithm [T| 


The notation and D 2 used in Algo, [^stands for the adjoint of Di and D 2 , respectively, 
while the notation prox stands for the proximity operator. The proximity operator associ¬ 
ated to a convex, lower semi-continuous function ip from T-l (where "H denotes a real Hilbert 
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Algorithm 1 Algorithm for solving (131 


Initialization 

Set (T > 0 and r G 


0 . 


i+IIEf=,jA0-AIP+3<^ 

Set hM 

^[0] -^^JxNixN^ ylo] G 

— =1 =2 


XN'2 


For ^ = 0,1.... 

Gradient descents steps 
_ r{D*y^^ + 

-2t (/iM (fc) _ u;[^l (j, fc) log2 A(j, 

Forj = ji,. ..,i 2 

■) = ■) - ■) 

-2^( ( Ej (ii log2 O', k) - (fe))log2A0, ^)) 

^ ' aGA- 


Proximity operator step based on dc-i 

Proximity operator step for TV on h 

= y[^ + - 4^^) 

(pr0X(;,y,)||.||^^^ (<7-1 (y[^(fc), (fc)))) 


Proximity operator step for dc 2 

i^+ll = f - (pro^,2>^cJo-^ql%,k)))_ 
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space) to ]—oo, +oo] is defined at the point u € "H as prox^(it) = argmin^jg-^^ |||u— i;||^+(^(d). 
The proximity operators involved in Algo, [^have a closed-form expression (cf. (4^): 


(Vue M^) prox >||,||2 j{u) = max(0,1 — 


A 


It. 


(7||U|| 2 


(16) 


Moreover, according to 50 Proposition 2.8]. if C denotes a non-empty closed convex subset 
of and if ry > 0, 


^ _l_ v{Pc{u)-u) if dr<(u) > n 

Vdcy-) jf dc{u)<7] 


(17) 


for every u € M'^. For our purpose, C models the hyperplane constraints Ci and C 2 and 
P^i Pc 2 have a closed-form expression given in 51 , consequently prox^^^^^ and prox,^^^ 
also have a closed-form expression. Note that the projection onto the convex set Pc is the 
proximity operator of the indicator function Lq of a non-empty closed convex set C C K 
(i.e., Lc{x) = 0 if a: S C and -|-oo otherwise). 

All example of the result provided by the proposed procedure is displayed in Fig. (f) 
for A = 16 and 7?i = 7/2 = 1000. It yields a piece-wise constant estimate that very well 
reflects the true regularity displayed in Fig. (b). Notably, the obtained histogram is 
pronouncedly spiked and can easily be thresholded in order to determine a labeling for the 
Q = 2 regions. For this example, the strategy clearly outperforms the more classical TV 


procedure of Section 3.2 


3.4 Direct estimation of local regularity labels 


The common feature of the 
viding denoised (piece-wise 
pointwise regularity is then 
Ill this section, we propose 
directly from the estimates 
estimates of the areas with 
and histogram thresholding 


approaches of Sections |3.2| and 3.3 is that they aim at first pro¬ 
constant) estimates of h. The labeling of regions with constant 
performed a posteriori by thresholding of the global histogram, 
a TV-based algorithm that addresses the partitioning problem 
h obtained using (j^ with a priori fixed weights w and yields 
constant regularity without recourse to intermediary denoising 
steps. 


Partitioning problem. Formally, the problem consists in identifying the areas 
(f2q)i<g<Q of a domain that are associated with different values {f^q)i<q<Q of h, 

(Vg e . .,0})(V^Co e nj /i(xo)-^g 

where fig = and (Vg 7^ p), fig fl flp = 0 (by convention, fig < Pg+i). Most methods 
for solving the partitioning problem are either based on the resolution of a non-convex 
criterion or require specific initialization [52|-[5^. Here, we adopt the minimal partitions 
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technique proposed in 1561, which is based on a convex relaxation of the Potts model and 


consequently enables convergence to a global minimizer. According to 57 , our partitioning 
problem can be written as 


minimize 


Q 

E 

q=l 


Q 


iq{h{x))dx+x Per(n^ 

9=1 


subj. to 


U?=, n, = 

(VgT^p), n9nf2p = 


(18) 


where Per(f2^) measures the perimeter of region tq{h{^) denotes the negative log- 
likelihood of the estimated local regularity associated with region Qq, and the constraints 
ensure a non-overlapping partition of (f^ 9 )i< 9 <Q- The parameter y > 0 models the roughness 
of the solution. 

Problem formulation. Let ^ e {0, I < q < Q, denote a set of Q partition 

matrices, i.e., 0 is -^i ^ N 2 matrix with all entries equal to one. The discrete 


analogue of (181 is the Potts model, which is known to be NP-hard to solve. A convex 


relaxation, involving the total variation, is given by (35l57 


Q 


32 


Q-l 


’ *^1 q=l k€K, j=3l 9=1 

subj. to {'ikeK] 1> ... >^^Q(fc)-0 (19) 

where the weights are a priori fixed, are independent of k. w{j^k) = and satisfy the 

constraints (j^. The regularization parameter A > 0 impacts the number of areas created 
for each single label. When A is small, several unconnected areas can occur for a single label 
while the solution favors dense regions when A is large. A bound on the error of the solution 


of the convex relaxation (19) is provided in |57| (for the special case of two classes the 
solution coincides with the global minimizer of the Ising problem (^), It results that the 


solutions of the minimization problem (19), denoted 6 S {0, g^j-g binary matrices 


that encode the partition matrices £2 such that 


dq_,{k) - 0,{k) = 


if Qq{k) = 1, 

otherwise. 


( 20 ) 


Algorithmic solution. The functions involved in (19) are convex, lower semi-continuous 
and proper, but the total variation penalty and the hard constraints are not smooth. The 
algorithm proposed in |57| considers the use of PDEs. In (^ , it is based on a Arrow- 
Hurwicz type primal-dual algorithm but requires inner iterations and upper boundedness 
of the primal energy in order to improve convergence speed. Here we employ a proximal 
algorithm in order to avoid inner iterations. Note that in |^, a proximal solution was 
proposed for a related minimization problem in the context of disparity estimation. For 
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the same reasons as those discussed in Section |3-3| we propose a solution based on the 
FBPD algorithm [^, specifically tailored to the problem (19l. The iterations are detailed 
in Algo. I^and involve the hyperplane constraints 


c, = {%_,^i^)\o,-Ak)>e,{k), k€K} 


for every q = 1..... Q with ^o(^) — 1 0Q{k) = 0, € /C. The projections onto Cq, 

denoted Pq {•), have closed-form expressions given in |51j. 


Algorithm 2 Algorithm for solving (191 


Initialization 

Set (T>najLdTG 

For 5 = l,...,Q-l, 

For g = 1,.... Q - 1. 2 ^ G 

=q =q -q 

For ^ = 0,1.... 

Gradient descents steps 
Forg ^ 1,... ,Q - 1 

- L) - "^ 4 ^^ 


xN. 


=9 


v.=9 


=9+1 =9 

-r{Dly 

For g = 1,..., Q odd 


=9 

m 

=9 


D2y 

=9 


2 ,W 


L ^= 9 - 


Proximity operator step based on TV and on even Cq 
Forg = T... - 1 

= yiM ^ o-Di(2£[^+^l - £M) 

g^.['l=ay*l+<jD2(2£['+il-£M) 

#1 =^m+^(29['+‘l-£l'l) 


=9 


=9 


=9 =9 


-o-prox„-i;^l|.|l2 j (cT 


For g = 1,..., Q even 


=9 =9 


-9-1 -g 

(T ’ (7 . 
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Negative log-likelihood. In the present work we focus on the Gaussian negative log- 
likelihood. i.e.. 


J2 

4 k) log2 X{i, k)) = 


k) log2 U, k) - Atq) 


3=Jl 


2a^ 


where fiq and denote the mean value and the variance in the region , respectively. 


=9 


The a priori choice of (Mg)i<g<(3 likely to strongly impact the estimates {0 )i<q<Q- 
We therefore propose to alternate the estimation of {9 )i<q<Q-i and {f^q)i<q<Q- The values 
{l-^q)i<g<Q first initialized equi-distantly between the minimum and the maximum values 
of h. Then Algo. j2 is rnii until convergence and the values (/^q)i<q<Q are re-estimated on the 
estimated areas (Q )i<q<Q- We iterate until stabilization of the estimates. The variances 


are fixed to = | here. 


4 Segmentation performance assessment 


Performance of the proposed procedures are qualitatively and quantitatively assessed using 

both synthetic scale-free textures (using realizations of 2D multifractional Brownian fields, 

described in Section 2.4 with prescribed piecewise constant local regularity values) and 

real-world textures, chosen in reference databases. Sample size is set to N x N = 512 x 

512 (hence, A^i = A ^2 = 256 due to the decimation operation in the DWT). Analysis is 

conducted using a standard 2D-DWT with orthonormal tensor product Daubechies mother 

wavelets with N^p = 2 vaiiLshing moments and 4 decomposition levels, (^ 1 ,^ 2 ) = (1)4). The 

labeling solutions obtained with the algorithms based on local smoothing, TV denoising, 

joint estimation of regularity and weights and direct local regularity labels are referred 
S '^TV TVW RJVlS 

to as n , , n and Q , respectively. For the algorithms involving a histogram 

thresholding step, thresholds are automatically set at the local minima between peaks of 
(smoothed) histograms (and, when less local minima than desired labels are detected, at 
the position of the largest peak). 


4.1 Quantitative performance assessment 

The proposed local regularity based labeling procedures are first compared using 10 indepen¬ 
dent realizations of mnltifractional Bro\^mian fields with Q = 2 areas of constant regularities 
hi and h 2 , respectively, given by the ellipse model shown in Fig.[^(b) (/i 2 corresponding to 
the inside of the ellipse). Performance are evaluated for a large range of values of the regu¬ 
larization parameter A (respectively, standard deviation a for the local Gaussian smoothing 
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Figure 2: Percentage of misclassified pixels as a function of penalty parameter A (re- 

-^S 

spectively, standard deviation a for ): median (solid red) and upper and lower quartile 

-"-TV —TVW -~RMS 

(dashed red) for 10 realizations for £2 , ^ , £2 (from left to right, respec¬ 

tively). Subfigure (a): (/ii,/i 2 ) = (0.5,0.7); Subfigure (b): {hi. h^} = (0.6,0.7). 



Figure 3: Re.sults obtained with the different proposed solutions compared to a basic smoo^i- 
ing of h and the state-of-the-art texture segmentation approaches proposed in 
Top row: (/ii,fi 2 ) = (0.5,0.7); bottom row: (/?i,fi 2 ) = (0.6,0.7). 
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(a) Data / (b) (c) 36 (d) M (e) g (f) g (i) g (j) 

Figure 4: Results of the labeling obtained with the different proposed solutions when Q = 3. 



Figure 5: Results obtained with the different proposed solutions (in red) compared to a 
basic smoothing (in blue) as a function of Ah = /i 2 ~ 


based solution). To assess performance, misclassified pixel rate is evaluated as follows: The 
area with the largest median of estimated local regularity values is associated with the 
area of the original mask with largest regularity, the area with the second largest median 
of estimated local regularity values with the area of the original mask with second largest 
regularity value, and so on. Achieved results are reported in Fig. (a) (for hi = 0.5 and 
/i 2 = 0.7) and (b) (hi = 0.6 and ^2 = 0.7), respectively. 


'-'TV 

The three TV-based strategies (D 


;-TVW ;-RMS, 

i I and i I 


clearly outperform the local 


smoothing based solution Q over a large range of A. The local smoothing procedure yields 
at best pixel misclassification rates of 12% (Fig. (a)) and 28% (Fig. (b)). In contrast, 
the best results obtained with the TV-based strategies drop down to less than 7% (Fig. 
(a)) and 12% (Fig. [^(b)) of misclassified pixels. 


-~TVW 

Among the TV-based algorithms, Q , relying on the joint estimation of regularity and 
weights, is the least sensitive to the precise selection of the regularization parameter A and 
consistently yields the best performance. Notably, it outperforms all other procedures when 
the difference in regularity |h 2 — hi| is small (Fig. (b)). For segmentation, the performance 
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— TV TVW —TV 

of the TV procedure (Q ) is similar to that of , yet Q is more sensitive to the precise 

tuning of A and yields large errors when A is chosen too small or too large. For these cases, 
— TV — TVW 

n detects only one area, while Q still segments the texture into two areas. Solution 

R,lvIS '^TV 

n also is more robust to the tuning of A, compared to , yet shows slightly decreased 

~ r -RMS “ 

misclassihcation rates. Furthermore, U has the practical advantage of being the only 
solution that does not require the practically cumbersome step of thresholding histograms 
(hence avoiding the empirical tuning of binning and smoothing parameters, for instance). 

To illustrate, results obtained by each of the proposed procedures for the value of A (or 
cr) leading to a minimal classification error (marked with symbols in Fig. are reported 
in Fig. [^when considering one randomly selected realization of multifractional Brownian 
fields with Q = 2 areas of constant regularity defined by {hi,h 2 ) = (0.5,0.7) (top) and 
(^ 1 )^ 2 ) = (0.6, 0.7) (bottom), respectively. Fig. (a) shows the analyzed textures, illus¬ 
trating that the two texture areas can not be distinguished visually. The local smooth- 

ing based labeling results Q are clearly the poorest, both for {hi. h 2 ) = (0.5, 0.7) and 

“ -TV —TVW —RMS 

(^ 1 )^ 2 ) = (0.6,0.7). In contrast, all three TV-based solutions Q , U and Q yield 

—TVW ~ ~ ~ 

satisfactory performance. Solution Q achieves the lowest misclassification error and is 

also visually the most convincing in terms of segmented regions, at the price though of the 

—TV 

largest computational cost. Solution Q shows only slightly larger misclassification rates 
sand may be preferred in certain applications for its significantly smaller computational 
cost. Solution Q shows larger misclassification rate when the difference in regularity 
decreases. 


Moreover, the above analysis is complemented by the study of a situation with oi Q = 3 
areas (constant local regularities (hi,fi 25 ^ 3 ) = (0.2, 0.4, 0.7)), with a more complex ge¬ 
ometry, including notably corners. The regularity mask, texture and labeling results are 
illustrated in Fig. The local smoothing solution Q. fails to distinguish the three areas, 

yielding several disconnected domains, while all proposed TV based approaches correctly 

is TV -TVW 

and satisfactorily detect the three distinct areas. Solutions Q and £2 show similar 

^ R,1VIS 

results and again have slightly smaller classification error rates than £2 . None of the 

proposed procedures recovers the two sharp corners of the original regularity mask. In 
-TV -TVW 

particular, £2 and £2 yield segmented areas with pronouncedly smooth borders. 

-TVW -TV 

In conclusion, the lowest misclassification rates are achieved by £2 , followed by £2 

and £2 , and worst results are obtained with £2 . The quality of the solutions obtained 

with the different strategies are related to their complexity, reflected by computation times. 

-s , , , -TV 

While £2 is obtained in less than 1 second, the fastest TV based solution is £2 , with a 

~ r • r -RMS 

computational time of about 10 seconds in our experiment, followed by £2 with a 1- 
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—TVW 

minute computational time. The solution achieving the best performance, Q , further 
requires aromid 3 minutes of computational time. 

In order to evaluate more accurately the behaviour of the proposed method as a func¬ 
tion of Ah = h 2 — hi, additional experiments are conducted in which /i 2 = hi + Ah for 
Ah = {0.1,0.2,0.3, 0.4}. Average (over 10 realizations) pixel misclassification rates for two 
different situations, hi = 0.2 and hi = 0.5, are reported in Fig (the reported results are 
obtained with values of A that lead to best performance for the configurations). The results 
based on TV are displayed in hot colors (orange/pink/red) while the results obtained with 
the basic smoothing is displayed in blue. As expected, the larger Ah, the better the per¬ 
formance for all methods. Moreover, the results lead to the conclusion that the TV based 
approaches have superior performance and that the method that estimates simultaneously 
the weights and the Holder exponent, i.e. Q , yields overall smallest misclassification 
rates. 

4.2 Comparisons with state-of-the-art segmentation procedures 
and application to real-world textures. 

Comparisons with state-of-the-art segmentation procedures. Comparisons 
against two texture segmentation procedures, chosen because considered state-of-the-art 
in the dedicated literature, are now discussed. The first approach |36| relies on a multiscale 
contour detection procedure using brightness, color and texture (using textons) information 
followed by the computation of an oriented watershed transform. The second approach 
relies on Gabor coefficients as features followed by a feature selection procedure relying on 
a matrix factorization step. B,esults are reported in Fig. and and unambiguously 

show that such approaches lead to much poorer segmentation performance as compared to 
the proposed TV-based procedures and appear ineffective for the segmentation of scale-free 
textures. 

Real-world textures. To assess the level of generality of the TV-based segmentation 
procedures, their performance are quantified and compared on sample textures chosen ran¬ 
domly from a large database considered as reference in the dedicated literature, the Berkeley 
Segmentation DatasetQ No ground truth segmentation is available for that database. Sam¬ 
ples are shown in Fig. [^and Fig. [^together with local regularity TV based estimation and 
segmentation procedure outcomes. For comparison, results achieved with the two state-of- 
the-art approaches in and are also displayed. For both examples, we provide two 
types of results: a 2-label segmentation result {Q = 2, Fig. and top rows) and the 
G-labels segmentation results where, for the state-of-the-art methods, Q had been tuned 
empirically in order to yield the visually most convincing solution, and where the solutions 

^ https: //www.eecs.berkeley.edu/Researcb/Projects/CS/vision/bsds 
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with Q = 9 with Q = 5 fl 


Figure 6: Experiments on real data: 
Database. 


-'TVW —RMS 

with Q = ion with Q = 9Q with <5 = 3 
image extracted from the Berkeley Segmentation 


0 


bottom rows). For both 


— TV —TVW . —TV —TVW FI 

n and Q are deduced from and (Fig- 6 and 

examples, we observe that the outcome of a 2-label segrnentationls overall more satisfactory 

for the proposed TV solution than for the state-of-the-art approaches |36| and [^, which 

tend to fail to extract meaningful shapes. When the number of labels Q is increased, the 

behaviour of the methods 1361 and [7| improves, yet the proposed TV methods, and in par- 
-TV —TVW *— ' r 

ticular h and h . also lead to mut’li more informative segmentation results. Finally, 
—— 1-1 

— TVW 

the results achieved with Q are observed to be visually more satisfactory than those 
of the two other TV-based approaches, notably for the second example. This is a very 
satisfactory outcome regarding the general level of applicability of the proposed TV-based 
segmentation approaches, as there is no reason a priori to believe that the samples in the 
Berkeley Segmentation Dataset are perfectly scale-free textures. Further randomly selected 
examples from that database are available upon request. 
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Figure 7: Experiments on real data: image extracted from the Berkeley Segmentation 
Database. 









5 Conclusions and perspectives 


We proposed, to the best of our knowledge, the first fully operational nonparametric tex¬ 
ture segmentation procedures that rely on the concept of local regularity. The segmentation 
procedures was designed for the class of piecewise constant local regularity images. The orig¬ 
inality of the proposed approach resided in the combination of wavelet leaders based local 
regularity estimation and proximal solutions for minimizing the convex criterion underlying 
the segmentation problem. Three original and distinct proximal solutions were proposed, 
all relying on a total variation penalization and proximal based resolution: TV denoising 
of local regularity estimates followed by thresholding for label determination; TV penal¬ 
ized joint estimation of local regularity and estimation weights followed by thresholding for 
label determination; direct labeling of local regularity estimates using a TV based partition¬ 
ing strategy. The performance of the procedures was validated using stochastic Gaussian 
model processes with prescribed region-wise constant local regularity and illustrated using 
realistic model images with real-world textures. All proposed labeling procedures yielded 
satisfactory results. They significantly improved over labeling based directly on (smoothed) 

local regularity estimates, at the price though of increased computational costs. Procedure 
-RMS 

h further avoided to devise a detailed procedure for histogram thresholding, yet yielding 
“ -TVW 

slightly poorer results compared to h 

When constant regularity areas were labeled, regularity can be re-estimated a posteriori 
by averaging within each area the X{a.k) prior to performing linear regressions [27| . 

The proposed procedures are currently being used for the analysis of biomedical textures, 
with encouraging preliminary results. Comparisons with alternative texture characterization 
features, such as local entropy rates, are under investigation. 

Future work will include investigating in how far the proposed approach can be adapted 
to handle different regularity models, such as images with piecewise smooth local regularity. 
Furthermore, the analysis of textures from real-world applications would benefit from sub¬ 
stituting piecewise constant Holder exponents with piecewise constant multifractal spectra, 
providing richer and more realistic models, at the price, yet, of more severe estimation and 
segmentation issues. 
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